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Abstract 

This  project  has  developed  and  analyzed  new  mathematical  algorithms  to  substantially  reduce 
the  complexity  of  simulating  and  optimizing  parametrically  dependent  systems  and  to  support 
decision-making  under  uncertainty.  Specifically,  this  research  has  advanced  the  state  of  the  art  in 
reduced  order  modeling  based  on  projections  and  on  the  discrete  empirical  interpolation  method 
(DEIM)  for  nonlinear  systems,  developed  new  adaptive  sampling  methods  for  optimization  of 
systems  with  uncertain  inputs,  devised  a  domain  decomposition  based  methods  to  systematically 
integrate  the  uncertainty  propagation  through  components  into  uncertainty  propagation  through  a 
systems  composed  of  these  components,  established  a  so-called  CUR  factorization  based  on  the 
DEIM  that  provides  a  low  rank  approximate  factorization  of  a  given  large  matrix  with  applications 
to  POD  model  reduction  and  analysis  of  data.  The  feasibility  of  our  algorithms  was  demonstrated 
on  a  number  of  problems  with  relevance  to  the  AirEorce. 


1  Introduction 


Effective  computational  tools  to  support  decision-making  under  uncertainty  have  become  essential 
in  the  design  and  operation  of  aerospace  systems.  Many  problems  of  relevance  to  the  Air  Eorce  are 
described  by  mathematical  models  of  high  complexity — involving  multiple  disciplines,  character¬ 
ized  by  a  large  number  of  parameters,  and  impacted  by  multiple  sources  of  uncertainty.  The  accu¬ 
rate  and  efficient  propagation  of  uncertainties  in  parameters  through  these  complex,  high  fidelity 
computational  models  is  a  significant  challenge.  This  project  has  developed  new  mathematical 
algorithms  and  analyses  to  substantially  reduce  the  complexity  of  simulating  and  optimizing  para¬ 
metrically  dependent  systems  and  of  propagating  uncertainty  through  systems.  In  particular,  this 
research  has  advanced  the  state  of  the  art  in  reduced  order  modeling  based  on  projections  and  on 
the  discrete  empirical  interpolation  method  (DEIM)  for  nonlinear  systems,  developed  new  adaptive 
sampling  methods  for  optimization  of  systems  with  uncertain  inputs,  devised  a  domain  decompo¬ 
sition  based  methods  to  systematically  integrate  the  uncertainty  propagation  through  components 
into  uncertainty  propagation  through  a  systems  composed  of  these  components,  established  a  so- 
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called  CUR  factorization  based  on  the  DEIM  that  provides  a  low  rank  approximate  factorization 
of  a  given  large  matrix  with  applications  to  POD  model  reduction  and  analysis  of  data.  The  fea¬ 
sibility  of  algorithms  was  demonstrated  on  a  number  of  problems  with  relevance  to  the  AirForce. 
This  report  highlights  several  of  our  research  results.  More  details  can  be  found  in  the  publications 
[1,  2,  3,  4,  5,  7,  8,  9,  10,  11,  13]  and  theses  [6,  12]  that  resulted  from  this  research. 

Sections  2.1  to  2.4  summarize  our  advancements  to  the  state  of  the  art  in  reduced  order  modeling 
based  on  projections  and  on  the  discrete  empirical  interpolation  method  (DEIM)  for  nonlinear  sys¬ 
tems.  Section  2. 1  is  concerned  with  the  application  of  the  Discrete  Empirical  Interpolation  Method 
(DEIM)  to  systems  obtained  by  finite  element  discretizations  of  partial  differential  equations,  and 
to  a  class  of  parameterized  system  that  arises,  e.g.,  in  shape  optimization.  Careful  analysis  and 
reorganizations  of  computations  lead  to  evaluations  of  the  DEIM  reduced  order  models  (ROMs) 
that  are  substantially  faster  than  those  of  the  standard  projection  based  ROMs.  Additional  gains 
are  obtained  with  the  DEIM  ROMs  when  one  has  to  compute  derivatives  of  the  model  with  respect 
to  the  parameter. 

Traditionally,  projection  based  reduced  order  models  use  only  one  subspace  for  all  parameters. 
To  obtain  an  accurate  ROM,  this  can  require  relatively  large  dimensional  subspaces,  and,  conse¬ 
quently,  result  in  a  less  efficient  ROM.  The  approach  in  Section  2.2  divides  the  parameter  space 
into  regions  and  for  each  region  develops  a  different  ROM.  The  regions  in  parameter  space  are  de¬ 
termined  using  machine  learning  methods  in  the  offline  computational  phase  via  clustering.  This 
approach  can  lead  to  significantly  smaller  and  therefore  computationally  cheaper  ROMs. 

Section  2.3  summarizes  our  development  of  an  efficient  approach  for  a-posteriori  error  estimation 
for  POD-DEIM  reduced  nonlinear  dynamical  systems.  Our  a-posteriori  error  estimator  can  be 
efficiently  decomposed  in  an  offline/online  fashion  and  is  obtained  by  a  one  dimensional  auxiliary 
ODE  during  reduced  simulations. 

The  application  of  the  DEIM  requires  that  a  component  of  the  system  nonlinearity  depends  only 
on  a  few  components  of  the  vector  of  arguments,  and  it  requires  explicit  knowledge  of  this  depen¬ 
dence.  If  DEIM  based  model  reduction  is  applied  to  existing  codes,  determining  this  dependency 
can  be  tedious  and  error  prone  task.  Section  2.4  describes  how  standard  techniques  of  automatic 
differentiation  (AD)  can  be  adapted  to  producing  such  a  code  and  hence  can  facilitate  the  use  of 
DEIM  by  non-experts. 

The  methods  in  sections  2.5  and  2.6  are  concerned  with  efficient  uncertainty  propagation  through 
coupled  systems,  and  efficient  optimization  of  PDF  systems  with  uncertain  parameters.  Thus, 
while  the  methods  summarized  in  sections  2.1  to  2.4  deal  with  complexity  reduction  in  simulation 
of  general  parameterized  systems,  the  approaches  in  sections  2.5  and  2.6  integrate  the  statistical 
distributions  of  the  uncertain  parameters. 

We  have  developed  a  rigorous  framework  based  on  trust-region  algorithms  and  adaptive  stochastic 
collocation  for  the  efficient  solution  of  optimization  problems  governed  by  partial  differential  equa¬ 
tions  with  uncertain  parameters.  Sparse  grid  collocation  or  other  sampling  methods  lead  to  inexact 
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objective  functions  and  gradients  in  the  optimization  problems.  Our  new  trust-region  framework 
rigorously  manages  the  use  of  this  inexact  problem  information.  As  a  result,  we  dramatically  re¬ 
duce  the  number  of  samples  required,  which  leads  to  a  dramatic  reduction  in  the  computational 
cost.  Moreover,  numerical  results  indicate  that  our  new  algorithm  rapidly  identifies  the  stochas¬ 
tic  variables  that  are  relevant  to  obtaining  an  accurate  optimal  solution.  When  the  number  of 
such  variables  is  independent  of  the  dimension  of  the  stochastic  space,  the  algorithm  exhibits  near 
dimension-independent  behavior. 

Decomposition  of  a  system  into  subsystems  or  component  parts  is  another  strategy  to  manage 
complexity.  We  have  developed  a  domain  decomposition  approach  for  uncertainty  analysis.  Our 
approach  sketched  in  Section  2.6  generates  computational  models  for  the  propagation  of  uncer¬ 
tainty  through  the  subsystems.  These  subsystem  models  are  computed  independently  in  an  “of¬ 
fline”  phase  and  they  account  for  propagation  of  uncertainly  through  interface  conditions,  which 
is  important  for  the  coupled  systems.  Once  the  subdomain  models  are  constructed,  importance 
sampling  is  used  to  couple  the  subdomain  models  and  analyze  the  propagation  of  uncertainty  in 
the  global  system.  This  last  step  is  performed  in  an“online”  phase  and  does  not  require  expensive 
simulations. 

Finally,  Section  2.7  summarizes  our  development  and  analyzes  of  a  new  CUR  matrix  factorization 
for  matrices  built  from  large-scale  data  sets.  The  CUR  offers  two  advantages  over  the  SVD:  when 
the  data  in  A  are  sparse,  so  too  are  the  C  and  R  matrices,  unlike  the  matrices  of  singular  vectors; 
and  the  columns  and  rows  that  comprise  C  and  R  are  representative  of  the  data  (e.g.,  sparse, 
nonnegative,  integer  valued,  etc.).  Our  algorithm  is  novel  in  that  it  only  requires  one  pass  through 
the  data  A  and  only  requires  a  few  columns  of  A  to  be  in  memory  at  any  given  step  of  the  algorithm. 
This  is  ideal  for  POD  based  model  reduction  schemes  since  the  high  dimensional  trajectory  need 
not  be  stored. 


2  Accomplishments 


This  rsection  highlights  several  of  our  research  results.  More  details  can  be  found  in  the  publica¬ 
tions  [1,  2,  3,  4,  5,  7,  8,  9,  10,  11,  13]  and  theses  [6,  12]  that  resulted  from  this  research. 


2.1  Application  of  the  Discrete  Empirical  Interpolation  Method  to  Reduced 
Order  Modeling  of  Nonlinear  and  Parametric  Systems 

Projection  based  methods  lead  to  reduced  order  models  (ROMs)  with  dramatically  reduced  num¬ 
bers  of  equations  and  unknowns.  However,  for  nonlinear  or  parametrically  varying  problems  the 
cost  of  evaluating  these  ROMs  still  depends  on  the  size  of  the  full  order  model  and  therefore  is 
still  expensive.  The  Discrete  Empirical  Interpolation  Method  (DEIM)  further  approximates  the 
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nonlinearity  in  the  projeetion  based  ROM.  The  resulting  DEIM  ROM  nonlinearity  depends  only 
on  a  few  eomponents  of  the  original  nonlinearity.  If  eaeh  eomponent  of  the  original  nonlinearity 
depends  only  on  a  few  eomponents  of  the  argument,  the  resulting  DEIM  ROM  ean  be  evaluated 
effieiently  at  a  eost  that  is  independent  of  the  size  of  the  original  problem.  Eor  systems  obtained 
from  finite  difference  approximations,  the  fth  component  of  the  original  nonlinearity  often  depends 
only  on  the  zth  component  of  the  argument.  This  is  different  for  systems  obtained  using  finite  ele¬ 
ment  methods,  where  the  dependence  is  determined  by  the  mesh  and  by  the  polynomial  degree  of 
the  finite  element  subspaces.  Our  paper  [AHS13]  describes  two  approaches  of  applying  DEIM  in 
the  finite  element  context,  one  applied  to  the  assembled  and  the  other  to  the  unassembled  form  of 
the  nonlinearity.  We  carefully  examine  how  the  DEIM  is  applied  in  each  case,  and  the  substantial 
efficiency  gains  obtained  by  the  DEIM.  In  addition,  we  demonstrate  how  to  apply  DEIM  to  obtain 
ROMs  for  a  class  of  parameterized  system  that  arises,  e.g.,  in  shape  optimization.  The  evaluations 
of  the  DEIM  ROMs  are  substantially  faster  than  those  of  the  standard  projection  based  ROMs. 
Additional  gains  are  obtained  with  the  DEIM  ROMs  when  one  has  to  compute  derivatives  of  the 
model  with  respect  to  the  parameter. 

The  following  example  is  a  semillinear  elliptic  reaction-advection-diffusion  partial  differential 
equation  (PDE)  posed  on  O  =  (0, 18)  x  (0,9)  x  (0,9).  A  concentration  of  w  =  0.2  is  imposed 
on  {0}  X  [3,6]  X  X  [3,6].  We  want  to  represent  the  concentration  (the  solution  of  the  semillinear 
elliptic  PDE)  for  varying  parameter  (ln(A) ,  E)  in  the  reaction  term.  A  finite  element  mesh  (Mesh  2) 
of  the  domain  Q  and  concentration  corresponding  to  parameter  (ln(A),E)  =  (6.4,0.11)  are  shown 
in  Eigure  1 . 


Eigure  1 :  Coarse  finite  element  mesh  of  the  domain  Q  (left  plot)  and  solution  of  the  semillinear 
reaction-advection-diffusion  PDE  with  parameter  (ln(A),E)  =  (6.4,0.11)  (right  plot). 
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piecewise  linear  FEM 

assembled  unassembled 


Figure  2:  The  tetrahedra  that  contain  DEIM  points  for  piecewise  linear  finite  elements  are  shown. 


Figure  2  shows  the  tetrahedra  in  Mesh  2  that  contain  a  node  corresponding  to  a  DEIM  point.  The 
plots  in  the  left  column  correspond  to  the  DEIM  applied  to  the  assembled  form  of  the  nonlinearity. 
The  expense  of  evaluating  the  DEIM  reduced  order  model  is  proportional  to  the  number  of  nodes 
the  DEIM  indices  are  connected  to.  Eigure  2  already  indicates  that  the  unassembled  form  of  DEIM 
is  more  efficient.  This  is  confirmed  by  the  timing  results  in  Table  1.  Table  1  clearly  show  the 
time  saving  of  DEIM-POD  receded  order  models  over  naive  POD  reduced  order  model  for  both 
piecewise  linear  {p  =  1)  and  piecewise  quadratic  {p  =  2)  finite  elements.  Using  the  unassembled 
form  of  DEIM  leads  to  additional  gains. 


Polynomial  degree 

P=  1 

P 

=  2 

Mesh  number 

I 

2 

3 

1 

2 

Pull 

1.78  (4) 

10.60  (3) 

43.30  (3) 

7.80  (3) 

185.00  (3) 

POD 

1.28  (4) 

8.04  (3) 

23.80  (3) 

4.12(3) 

38.80  (3) 

POD-DEIM 

0.15  (4) 

0.10(3) 

0.21  (4) 

0.21  (3) 

0.40  (3) 

POD-DEIM-u 

0.16(9) 

0.07  (4) 

0.10(4) 

0.01  (4) 

0.18(4) 

Table  1:  The  computing  times  (in  sec)  and  the  number  of  Newton  iterations  (in  parenthesis)  needed 
to  solve  the  full  order  model,  the  POD  reduced  order  model,  the  POD-DEIM  reduced  order  model, 
and  the  POD-DEIM-u  (unassembled)  reduced  order  model  for  different  grid  levels  and  linear  {p  = 
1)  and  quadratic  (p  =  2)  finite  elements. 


2.2  Localized  DEIM 

We  have  proposed  the  Eocalized  DEIM  (EDEIM)  as  way  to  improve  the  efficiency  of  DEIM  ROMs 
for  nonlinear  parameterized  PDEs  [PBWB13].  Whereas  DEIM  projects  the  nonlinear  term  onto 
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one  global  subspace,  our  LDEIM  method  computes  k  local  subspaces,  each  tailored  to  a  particu¬ 
lar  region  of  characteristic  system  behavior.  Then,  depending  on  the  current  state  of  the  system, 
LDEIM  selects  an  appropriate  local  subspace  for  the  approximation  of  the  nonlinear  term.  In  this 
way,  the  dimensions  of  the  local  DEIM  subspaces,  and  thus  the  computational  costs,  remain  low 
even  though  the  system  might  exhibit  a  wide  range  of  behaviors  as  it  passes  through  different 
regimes.  LDEIM  uses  machine  learning  methods  in  the  offline  computational  phase  to  discover 
these  regions  via  clustering.  In  our  implementation,  the  snapshots  are  clustered  with  A:-means.  Lo¬ 
cal  DEIM  approximations  are  then  computed  for  each  cluster.  In  the  online  computational  phase, 
machine-learning-based  classification  procedures  select  one  of  these  local  subspaces  adaptively  as 
the  computation  proceeds.  We  use  a  nearest  neighbor  classifier.  The  classification  can  be  achieved 
using  either  the  system  parameters  or  a  low-dimensional  representation  of  the  current  state  of  the 
system  obtained  via  feature  extraction. 

The  LDEIM  approach  is  demonstrated  for  a  reacting  flow  example  of  a  premixed  H2-Air  flame. 
The  physics  are  modeled  with  the  one-step  reaction  mechanism 

2H2  +  02^2H20, 

where  H2  is  the  fuel,  O2  is  the  oxidizer,  and  H2O  is  the  product.  The  evolution  of  the  flame  in  the 
domain  Q.  is  given  by  the  nonlinear  advection-diffusion-reaction  equation 

KAy-wVy-l-f(y,A/)  =0,  (1) 

where  y  =  [yYi2,yo2^y^20^TY  contains  the  mass  fractions  of  species  H2,02,  and  H2O  and  the 
temperature.  The  constants  k  =  2.0cm^/sec  and  w  =  50cm/sec  are  the  molecular  diffusivity  and 
the  velocity  of  the  velocity  field  in  xi  direction,  respectively.  The  geometry  of  the  problem  is  shown 
in  Eigure  3.  We  prescribe  homogeneous  Dirichlet  boundary  conditions  on  the  mass  fractions  on  E 1 
and  r3,  and  homogeneous  Neumann  conditions  on  temperature  and  mass  fractions  on  r4,r5,  and 
r6.  We  have  Dirichlet  conditions  on  r2  with  =  0.0282,  yo2  =  0.2259,  yH20  =  0,yr  =  950K, 

and  on  Fi ,  r3  with  yr  =  300K. 


Xi 


Eigure  3:  The  spatial  domain  of  the  reacting  flow  simulation. 


The  nonlinear  reaction  source  term  f(y,A/  =  [/H2(y,A'),/o2(y,A'),/H2o(y,A'),/r(y,A')]^  in  (1)  has 
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the  components 


/i(y,A')  =  -Vi(riH2yH2)^(T102yo2)A'iexp  ,  /  =  H2,02,H20 

/r(y,)u)  =  2/H2o(y,)u), 

where  V;  and  T],-  are  constant  parameters,  R  =  8.314472J/(mol  K)  is  the  universal  gas  constant,  and 
Q  =  9800K  is  the  heat  of  reaction.  The  parameters  /j  =  (^ui  ,/r2)  £  ®  with  (D  =  [5.5e+l  1 , 1 .5e+03]  x 
[1.5e+13,9.5e+03]  are  the  pre-exponential  factor  and  the  activation  energy,  respectively.  The  equa¬ 
tion  (1)  is  discretized  using  the  finite  difference  method  on  a  73  x  37  grid  leading  to  =  10, 804 
degrees  of  freedom.  The  resulting  nonlinear  system  of  equations  is  solved  with  the  Newton  method. 
The  snapshots  are  computed  for  the  parameters  on  a  50  x  50  equidistant  grid  in 

Figure  4  compares  DEIM,  parameter-based  LDEIM  with  splitting  and  clustering,  and  state-based 
EDEIM.  Eor  splitting,  we  set  the  tolerance  £  to  le-08,  which  is  about  two  orders  below  what  DEIM 
achieves.  Eor  the  parameter-based  EDEIM  with  clustering,  the  number  of  clusters  is  set  to  5.  More 
clusters  lead  to  an  unstable  behavior.  Eor  the  state-based  EDEIM,  however,  we  set  the  number  of 
clusters  to  15.  The  state-based  EDEIM  uses  a  point-based  feature  extraction  with  M  =  5  dimen¬ 
sions.  In  all  cases,  we  have  40  POD  modes.  In  Eigure  4,  we  see  that  the  results  of  the  parameter- 
based  EDEIM  with  clustering  do  not  improve  after  about  10  DEIM  modes.  Again,  the  clustering 
becomes  unstable.  However,  this  is  not  the  case  for  state-based  EDEIM  which  achieves  about  two 
orders  of  magnitude  better  accuracy  than  DEIM.  The  same  holds  for  the  splitting  approach. 


Eigure  4:  Eor  the  reacting  flow  example,  the  comparison  between  DEIM,  parameter-based  EDEIM, 
and  state-based  EDEIM.  The  results  are  shown  for  40  POD  and  20  DEIM  modes.  State-based 
EDEIM  achieves  about  two  orders  of  magnitude  better  accuracy  than  DEIM  for  the  same  number 
of  DEIM  modes. 
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2.3  A-Posteriori  Error  Estimation  for  DEIM 


We  have  introdueed  an  effieient  approaeh  for  a-posteriori  error  estimation  for  POD-DEIM  redueed 
nonlinear  dynamieal  systems  [WSH12].  The  eonsidered  nonlinear  systems  may  also  inelude  time 
and  parameter-affine  linear  terms  as  well  as  parametrieally  dependent  inputs  and  outputs.  The  re- 
duetion  proeess  involves  a  Galerkin  projeetion  of  the  full  system  and  approximation  of  the  system’s 
nonlinearity  by  the  DEIM  method  [Chaturantabut  &  Sorensen  (2010)].  The  proposed  a-posteriori 
error  estimator  can  be  efficiently  decomposed  in  an  offline/online  fashion  and  is  obtained  by  a  one 
dimensional  auxiliary  ODE  during  reduced  simulations.  Key  elements  for  efficient  online  compu¬ 
tation  are  partial  similarity  transformations  and  matrix  DEIM  approximations  of  the  nonlinearity 
Jacobians.  The  theoretical  results  are  illustrated  by  application  to  an  unsteady  Burgers  equation 
and  a  cell  apoptosis  model. 

The  key  ideas  of  our  error  estimation  procedure  is  an  application  of  a  Gronwall-like  differential 
inequality  for  the  error  norm  ||y(t)  —  yr(t)  ||  (full  vs.  reduced  state) ,  estimation  of  local  logarithmic 
Eipschitz  constants  and  an  estimation  of  the  DEIM  approximation  error  for  a  ROM  of  order  m 
using  a  higher  order  DEIM  approximation  of  order  m  +  m' .  Essentially,  we  treat  the  higher  order 
approximation  heuristically  as  representing  the  true  state.  We  use  these  tools  to  provide  a-posteriori 
error  bounds  which  can  be  computed  along  with  the  reduced  system  trajectory  at  negligible  cost. 
Only  a  simple  ID  ODE  must  be  propagated  along  with  the  reduced  state  calculation. 


Unf-L2  absolute  e 


order 


Estimated  absolute  error 


Eigure  5:  (Eeft)  Absolute  error  of  DEIM  Approximation  (Right)  Demonstration  of  Error  Estimate 
for  varying  m  and  m' 


The  graphs  in  Eigure  6  illustrate  the  effectiveness  of  our  error  estimator.  The  graph  on  the  left 
shows  the  absolute  error  of  the  reduced  state  from  the  DEIM  reduced  model  with  the  true  state. 
The  right  graph  is  a  surface  plot  showing  the  error  estimate  using  the  m  +  m'  model  in  place  of  the 
true  state.  The  axis  labeled  “order”  shows  the  error  as  the  model  order  m  is  increased.  The  axis 
labeled  “error  order”  shows  the  estimated  error  as  m'  increases.  The  fact  that  the  surface  quickly 
flattens  out  in  the  m'  direction  verifies  that  m'  does  not  need  to  be  large,  it  typically  adequate  for 
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accurate  error  estimation  at  values  between  5  and  10.  This  plot  is  for  a  ID  Burgers’  equation  but 
we  obtained  virtually  identical  behavior  for  a  60,000  variable  2D  Reaction-Diffusion  model  Cell 
Apoptosis. 


2.4  Automating  DEIM 

As  mentioned  previously,  the  Discrete  Empirical  Interpolation  Method  (DEIM)  is  based  upon  a 
modification  to  proper  orthogonal  decomposition  which  is  designed  to  reduce  the  computational 
complexity  for  evaluating  the  reduced  order  nonlinear  term.  The  DEIM  approach  is  based  upon  an 
interpolatory  projection  and  only  requires  evaluation  of  a  few  selected  components  of  the  original 
nonlinear  term.  Thus,  implementation  of  the  reduced  order  nonlinear  term  requires  a  new  code  to 
be  derived  from  the  original  code  for  evaluating  the  nonlinearity.  We  have  developed  a  methodol¬ 
ogy  for  automatically  deriving  a  code  for  the  reduced  order  nonlinearity  directly  from  the  original 
nonlinear  code.  The  methodology  is  an  adaptation  of  standard  techniques  of  automatic  differenti¬ 
ation  (AD).  This  capability  removes  the  possibly  tedious  and  error  prone  task  of  producing  such  a 
code  by  hand  and  hence  can  facilitate  the  use  of  DEIM  by  non-experts. 

We  have  a  working  prototype  system  that  establishes  proof  of  concept  [CS12].  We  were  able 
to  automatically  produce  a  reduced  order  code  from  the  original  code  for  2-phase  miscible  flow 
simulation.  This  automatically  derived  code  was  able  to  completely  reproduce  the  simulation 
results  from  the  hand  coded  ROM  with  essentially  no  degradation  in  computation  time.  The  idea 
is  to  first  obtain  a  computational  graph  representation  of  the  original  code  and  then  this  graph  is 
“pruned”  to  retain  only  the  graph  necessary  to  compute  the  code  to  evaluate  the  DEIM  selected 
nodes  (which  correspond  to  the  DEIM  indices  of  the  desired  components  of  the  nonlinearity).  The 
pruned  graph  is  then  used  in  a  reverse  process  to  generate  the  code  for  the  ROM.  Using  standard 
AD  techniques,  one  could  also  generate  the  Jacobian  of  the  ROM  automatically. 


2.5  Trust-Region  Algorithms  with  Adaptive  Stochastic  Collocation  for  PDE 
Optimization  under  Uncertainty. 

We  have  improved  our  trust-region  adaptive  sparse  grid  framework  for  the  solution  of  optimization 
problems  governed  by  partial  differential  equations  (PDEs)  with  random  coefficients.  While  previ¬ 
ously  adaptive  sparse  grids  were  only  used  to  generate  the  models  used  to  commute  the  step  inside 
the  trust-region  framework  [7],  we  now  also  allow  the  approximation  of  function  values  using 
adaptive  sparse  grids  [8].  This  results  in  substantial  computational  savings.  With  our  algorithmic 
improvements,  numerical  solutions  that  previously  required  high  performance  parallel  computers^ 
can  now  be  run  on  a  desktop  [8]. 

'The  computations  in  [7]  were  carried  out  on  RedSky,  an  institutional  computing  cluster  at  Sandia  National  Labo¬ 
ratories.  The  cluster  is  built  on  a  3D  toroidal  QDR  InfiniBand  interconnect  and  provides  2,816  compute  nodes,  each 
with  8  cores,  2.93  GHz  Nehalem  X5570  processors,  and  12  GB  RAM  per  compute  node. 
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The  numerical  solution  of  optimization  problems  governed  by  PDEs  with  random  coefficients  po¬ 
tentially  requires  a  the  large  number  of  deterministic  PDE  solves  at  each  optimization  iteration. 
Our  algorithm  for  solving  such  problems  based  on  a  combination  of  adaptive  sparse-grid  collo¬ 
cation  for  the  discretization  of  the  PDE  in  the  stochastic  space  and  a  trust-region  framework  for 
optimization  and  fidelity  management  of  the  stochastic  discretization.  The  overall  algorithm  adapts 
the  collocation  points  based  on  the  progress  of  the  optimization  algorithm  and  the  impact  of  the 
random  variables  on  the  solution  of  the  optimization  problem.  It  frequently  uses  few  collocation 
points  initially  and  increases  the  number  of  collocation  points  only  as  necessary,  thereby  keeping 
the  number  of  deterministic  PDE  solves  low  while  guaranteeing  convergence.  Our  trust-region 
framework  now  allows  adaptation  to  build  the  model,  as  well  as  to  approximate  the  objective  func¬ 
tion  [8].  Our  trust-region  framework  can  be  used  in  a  wide  range  of  applications,  and  extends  the 
state-of-the-art  in  trust  region  methods  for  problems  with  inexact  function  and  derivative  infor¬ 
mation.  Currently  an  error  indicator  is  used  to  estimate  gradient  errors  due  to  adaptive  stochastic 
collocation,  but  other  adaptive  discretizations  could  be  used  as  well. 

The  following  example  is  motivated  by  the  optimal  control  of  direct  field  acoustic  testing  (DEAT). 
An  important  objective  in  DEAT  is  to  accurately  shape  sound  pressure  fields  in  a  region  of  interest 
by  using  acoustic  sources,  such  as  loudspeakers,  located  away  from  and  directed  at  the  region  of 
interest.  When  the  refractive  index  of  the  medium  in  the  region  of  interest  is  random,  we  may 
assume  that  the  governing  physics  are  modeled  by  the  stochastic  Helmholtz  equation 

—Au{y,x)  —  k^{\+<5t{y,x))^u{y,x)—z{x)  \/{y,x)&TxD  (2) 

with  Robin  boundary  conditions  ^u{y,x)  =  iku{y,x)  for  all  {y,x)  G  E  x  dD.  Here  x  is  the  spatial 
variable  in  the  spatial  domain  D  and  y  is  the  random  variable.  We  consider  an  idealized  DEAT 
example  in  two  spatial  dimensions,  where  the  domain  is  D  =  (—5,5)^.  The  goal  of  this  control 
problem  is  to  match  the  wave  pressure  m  to  a  desired  wave  pressure  w  G  L?'{D;C)  in  the  disk 
Dr  =  {x  G  D  :  ||jc||2  <  R}  C  D.  We  apply  distributed  controls  in  Dc,  which  is  the  union  of  the 
regions  occupied  by  20  loudspeakers  arranged  in  a  circle  exterior  to  Dr.  Eigure  6  depicts  the 
experimental  setup. 
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Figure  6:  Left:  A  sketch  of  the  computational  domain  with  the  region  of  interest  Dr  and  20  loud¬ 
speakers  spaced  uniformly  along  a  circle  surrounding  Dr.  The  red  (lighter)  color  denotes  uncer¬ 
tainty  in  the  refraction  index  inside  Dr  and  uncertainty  in  the  loudspeaker  enclosures  (including 
thickness  and  sound  speed).  The  blue  (darker)  color  stands  for  the  loudspeaker  control  regions, 
whose  union  is  denoted  by  Dc-  Right:  The  real  component  of  the  desired  state,  w. 


In  this  example  we  study  three  sources  of  uncertainty:  refractive  index,  loudspeaker  thickness, 
and  speed  of  sound  in  each  loudspeaker.  We  assume  that  the  stochastic  refractive  index  is  of  the 
form  £(y,x)  =  Em=i  ^m{x)ym  and  that  its  covariance  is  given  by  a  Matern  covariance  functions.  The 
series  representation  of  £  is  motivated  by  a  truncated  Karhunen-Loeve  (KL)  expansion.  In  total,  we 
have  dim  =  40  +  M  random  variables,  40  random  loudspeaker  parameters  (enclosure  thicknesses 
and  wave  number  for  each  speaker)  and  M  random  variables  representing  the  stochastic  refractive 
index  £.  For  additional  details  see  Section  5.2  in  [KHRW14].  Figure  7  shows  the  optimal  solution 
for  {dim  =  60  random  variables  (i.e.,  M  =  20). 


Figure  7:  Left:  The  real  parts  of  the  computed  optimal  controls.  Center:  The  real  part  of  the 
expected  value  of  the  optimal  state,  restricted  to  the  region  of  interest.  Right:  The  real  part  of  the 
standard  deviation  of  the  optimal  state. 
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Table  2  gives  the  computational  cost  of  our  algorithm  as  the  stochastic  dimension  increases  from 
dim  =  42  to  dim  =  80  and  it  illustrates  the  power  of  our  algorithm.  Even  for  the  largest  dimension, 
dim  =  80,  the  problem  can  now  be  solved  using  a  desktop.  In  contrast,  the  40-dimensional  problem 
using  the  previous  version  of  our  algorithm  [KHRW13]  required  1,804,001  collocation  points 
for  the  evaluation  of  the  high-fidelity  objective  function,  while  the  80-dimensional  problem  now 
requires  only  311  points  due  to  the  adaptive  objective  function  evaluations.  We  also  note  that  in  this 
example  many  of  the  random  variables  are  not  important  for  the  optimization  context.  However, 
it  is  not  clear  which  one  are  the  important  ones.  Our  algorithm  is  able  to  automatically  detect  the 
important  random  variables  and  automatically  zooms  in  on  approximately  10  (out  of  80)  stochastic 
dimensions  that  are  relevant  to  achieving  objective  function  and  gradient  consistency  conditions 
required  by  the  trust-region  framework.  As  a  result  it  vastly  reduces  the  effective  problem  size 
and,  in  this  example,  shows  a  convergence  behavior  that  is  nearly  independent  of  the  stochastic 
dimension. 


dim 

PDE  Solves 

CPobj 

CPgrad 

Obj.  Value 

42 

11,543 

145 

145 

5.2542 

44 

32,739 

233 

481 

5.2637 

46 

60,617 

243 

1,453 

5.2641 

48 

79,221 

247 

2,961 

5.2641 

50 

90,157 

251 

4,569 

5.2641 

60 

100,911 

271 

7,621 

5.2641 

70 

103,979 

291 

8,233 

5.2641 

80 

105,607 

311 

8,253 

5.2641 

Table  2:  Computational  cost  of  the  fully  adaptive  trust-region  algorithm  applied  to  the  Helmholtz 
control  example.  Here  dim  is  the  stochastic  dimension,  PDE  Solves  is  the  total  number  of  forward 
and  adjoint  PDE  solves,  CPobj  is  the  final  number  of  collocation  points  used  by  the  adaptive 
objective  function  scheme,  CPgrad  is  the  final  number  of  collocation  points  used  by  the  adaptive 
subproblem  (gradient)  model,  and  Obj.  Value  is  the  computed  value  of  the  objective  function  at 
termination. 


2.6  A  Domain-Decomposition  Approach  for  Uncertainty  Analysis. 

We  have  proposed  a  new  decomposition  approach  for  uncertainty  analysis  of  systems  governed 
by  PDEs  [EW14].  The  system  is  split  into  local  components  using  domain  decomposition.  Our 
domain-decomposed  uncertainty  quantification  (DDUQ)  approach  performs  uncertainty  analysis 
independently  on  each  local  component  in  an  “offline”  phase,  and  then  assembles  global  uncer¬ 
tainty  analysis  results  using  pre-computed  local  information  in  an  “online”  phase.  At  the  heart  of 
the  DDUQ  approach  is  importance  sampling,  which  weights  the  pre-computed  local  PDE  solutions 
appropriately  so  as  to  satisfy  the  domain  decomposition  coupling  conditions.  To  avoid  global  PDE 
solves  in  the  online  phase,  a  proper  orthogonal  decomposition  reduced  model  provides  an  efficient 
approximate  representation  of  the  coupling  functions. 
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The  DDUQ  approach  comprises  the  following  three  steps:  (1)  pre-step:  generate  POD  basis  for 
interface  functions;  (2)  offline  step:  generate  local  solution  samples  and  construct  surrogates  for 
the  coupling  functions;  (3)  online  step:  estimate  target  input  PDFs  using  surrogates  and  re-weight 
output  samples.  The  pre-step  is  cheap,  since  although  it  performs  standard  domain  decomposition 
iterations,  the  number  of  input  samples  is  small.  In  the  offline  stage,  expensive  computational  tasks 
are  fully  decomposed  into  subdomains,  i.e.,  there  is  no  data  exchange  between  subdomains.  The 
online  stage  is  relatively  cheap,  since  no  PDF  solve  is  required.  This  approach  is  summarized  in 
Figure  8,  where  communication  refers  to  data  exchange  between  subdomains. 


Pre-step  (cheap) 

•  Generate  interface  POD  basis. 

(Solve  PDEs,  has  communication,  but  a  few  samples) 


Offline  step  for  Di  (expensive) 

•  Generate  local  solution  samples. 

•  Construct  coupling  surrogates. 
(Solve  PDEs,  no  communication) 


Offline  step  for  Dm  (expensive) 

•  Generate  local  solution  samples. 

•  Construct  coupling  surrogates. 
(Solve  PDEs,  no  communication) 


Online  step  (cheap) 

•  Estimate  target  input  PDEs  using  coupling  surrogates. 

•  Re-weight  offline  output  samples. 

(No  PDE  solve,  has  communication) 


Figure  8:  Domain-decomposed  uncertainty  quantification  (DDUQ)  summary. 


In  the  offline  stage  of  DDUQ,  we  first  specify  a  proposal  input  PDF  from  which  the  samples  will 
be  drawn.  We  denote  the  proposal  input  PDF  by  where  is  the  vector  of  uncertain 

parameters  corresponding  to  subdomain  i  and  x,-  is  the  vector  of  interface  parameters  for  subdomain 
i.  The  proposal  input  PDF  must  be  chosen  so  that  its  support  is  large  enough  to  cover  the  support 
of  the  true  (unknown)  interface  parameters.  Provided  this  condition  on  the  support  is  met,  any 
proposal  input  PDF  can  be  used;  however,  the  better  the  proposal  (in  the  sense  that  it  generates 
sufficient  samples  over  the  support  of  the  target  input  PDF)  the  better  the  performance  of  the 
importance  sampling.  A  poor  choice  of  proposal  input  PDF  will  lead  to  many  wasted  samples  (i.e., 
requiring  many  local  PDE  solves). 

The  next  step,  still  in  the  offline  phase,  is  to  perform  uncertainty  quantification  on  each  local 
subdomain  i  independently.  This  involves  generating  a  large  number  A,  of  samples 
of  t;.(^/,x,),  where  the  superscript  (5)  denotes  the  5th  sample,  and  computing  the  local  solutions 

{m(v:,^-''\x-'*^)}  by  solving  the  deterministic  problem  for  each  sample.  Once  the  local  solutions  are 
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computed,  we  evaluate  the  local  outputs  of  interest  and  store  them. 

In  the  online  stage  of  DDUQ,  we  first  generate  Nq^  samples,  of  the  joint  PDF  of  inputs 

For  each  sample  we  use  the  domain  decomposition  iteration  to  evaluate  the  corre¬ 
sponding  interface  parameters.  After  the  above  process,  we  have  obtained  a  set  of  samples  drawn 
from  each  local  target  input  PDF  Ti:^.  x,).  We  next  estimate  each  local  target  input  PDF  from 

the  samples  using  a  density  estimation  technique,  such  as  kernel  density  estimation  (KDE). 

The  final  step  of  the  DDUQ  online  stage  is  to  use  importance  sampling  to  re- weight  the  outputs 
|y/  (^u  j  j  I  that  we  generated  by  PDF  solves  in  the  offline  stage.  The  weights,  wf\ 

are  computed  by  taking  the  ratio  of  the  estimated  target  PDF  to  the  proposal  PDF  for  each  local 
subdomain: 


z  =  1,...,M. 


We  can  show  that  the  probability  computed  from  these  weighted  samples 

j  j  I  ^  is  consistent  with  the  actual  distributions  of  the  outputs,  so 


that  under  certain  conditions. 


We  illustrate  the  DDUQ  approach  using  a  diffusion  problem  posed  on  a  2D  spatial  domain  with  two 
subdomains  with  randomness  in  the  permeability  coefficient,  the  Robin  coefficient,  and  the  source 
function.  The  PDFs  of  the  outputs  of  this  problem  are  shown  in  Figure  9,  where  we  see  that  as 
the  number  of  offline  samples  Aoff  increases  the  PDFs  generated  by  DDUQ  approach  the  reference 
PDFs  (generated  using  the  system-level  Monte  Carlo  simulation  with  Aref  =  10^  samples). 


Figure  9:  PDFs  of  the  outputs  of  interest  for  a  four-parameter  test  problem  using  the  diffusion 
equation. 
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2.7  CUR  Factorization 


We  have  developed  a  new  CUR  matrix  faetorization  based  upon  the  Diserete  Empirieal  Interpola¬ 
tion  Method  (DEIM)  [SE14].  A  CUR  faetorization  provides  a  low  rank  approximate  faetorization 
of  a  given  matrix  A  G  that  is  of  the  form  A  ^  CUR  where  C  =  A(:,  q)  G  is  a  subset  of 
the  eolumns  of  A  and  R  =  A(p, :)  G  is  a  subset  of  the  rows  of  A.  The  kxk  matrix  U  is  eon- 
strueted  to  assure  that  CUR  is  a  good  approximation  to  A.  Assuming  the  best  rank-k  singular  value 
deeomposition  (SVD)  A  ^  VSW^  is  available,  the  algorithm  uses  the  DEIM  points  q  =  DEIM(V) 
and  p  =  DEIM(W)  to  seleet  the  matriees  C  and  R. 

This  approximate  faetorization  is  nearly  as  aeeurate  as  the  best  rank  k  SVD  approximation  in  the 
sense  that 

||A-CUR||2=0(o,+i) 

where  O/t+i  is  the  first  negleeted  singular  value  of  A. 

The  CUR  faetorization  is  an  important  tool  for  matriees  built  from  large-seale  data  sets,  in  whieh 
ease  it  offers  two  advantages  over  the  SVD:  when  the  data  in  A  are  sparse,  so  too  are  the  C  and 
R  matriees,  unlike  the  matriees  of  singular  veetors;  and  the  eolumns  and  rows  that  eomprise  C 
and  R  are  representative  of  the  data  (e.g.,  sparse,  nonnegative,  integer  valued,  ete.).  The  following 
simple  example,  adapted  from  Mahoney  and  Drineas^,  illustrates  the  latter  advantage.  Suppose 
that  A  G  has  eolumns  eaeh  of  whieh  are  one  of  the  forms 


.ri 

,  or 

■  -1 

1  ■ 

.ri 

.  -^2  . 

2 

1 

1 

.  -^2  _ 

where  in  both  oases  xi  ~  A^(0, 1)  and  X2  ~  V(0,4^)  are  independent  normal  random  samples.  Thus 
the  eolumns  of  A  are  drawn  from  two  different  multivariate  normal  distributions.  Eigure  10  shows 
that  the  two  left  singular  veetors,  though  orthogonal  by  oonstruotion,  fail  to  represent  the  true 
nature  of  the  data;  in  oontrast,  the  first  two  eolumns  seleoted  by  the  DEIM-CUR  prooedure  give 
a  mueh  better  overall  impression  of  the  data.  While  trivial  in  this  two-dimensional  ease,  one  ean 
imagine  the  utility  of  sueh  approximations  for  high-dimensional  data. 

^CUR  matrix  decompositions  for  improved  data  analysis.  Proc.  Nat.  Acad.  Sci.,  106:697-702,  2009 
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Figure  10:  Comparison  of  singular  vectors  (left,  scaled,  in  red)  and  DEIM-CUR  columns  (right, 
in  blue)  for  a  data  set  drawn  from  two  multivariate  normal  distributions  having  different  principal 
axes. 


CUR-type  factorizations  had  their  origin  in  “pseudoskeleton”  approximations^  and  pivoted,  trun¬ 
cated  QR  decompositions^.  Numerous  prominent  recent  algorithms  instead  use  leverage  scores^. 
These  algorithms  all  rely  upon  initially  having  a  low  rank  SVD  approximation.  The  error  analysis 
of  these  algorithms  is  probabilistic  in  nature. 

Our  algorithm,  based  upon  DEIM,  is  entirely  deterministic  and  involves  few  (if  any)  parameters. 
We  have  developed  an  error  analysis  that  applies  to  a  broad  class  of  CUR  factorizations.  We  have 
also  proposed  a  novel  incremental  QR  algorithm  for  approximating  the  SVD.  This  approximate 
QR  algorithm  can  also  be  applied  to  efficiently  compute  leverage  scores  if  desired.  The  excellent 
performance  of  our  new  CUR  factorization  is  illustrated  on  several  examples.  These  results  indi¬ 
cate  that  the  DEIM-CUR  approach  generally  provides  far  superior  low  rank  approximations  to  the 
given  data  when  compared  to  existing  leverage  score  based  methods. 

The  incremental  low  rank  QR  factorization  has  important  consequences  for  POD  based  model 
reduction.  Here,  the  given  data  A  VT  where  the  columns  of  V  are  mutually  orthonormal  and 
T  is  rectangular.  Given  a  user  specified  tolerance  tol,  the  procedure  will  automatically  produce 
a  rank  k  factorization  that  achieves  ||A  —  VT||  <  to/||T||  where  V  is  order  nxk  while  T  is  order 

^S.  A.  Goreinov,  E.  E.  Tyrtyshnikov,  and  N.  L.  Zamarashkin.  A  theory  of  pseudoskeleton  approximations.  Linear 
Algebra  AppL,  261:1-21,  1997. 

W.  Stewart.  Eour  algorithms  for  the  efficient  computation  of  truncated  QR  approximations  to  a  sparse  matrix. 
Numer.  Math.,  83:313-323,  1999. 

^C.  Boutsidis  and  D.  R  Woodruff.  Optimal  CUR  matrix  decompositions,  arXiv.cs.DS:  1405.7910,  2014.  R  Drineas, 
M.  W.  Mahoney,  and  S.  Muthukrishnan.  Relative-error  CUR  matrix  decompositions.  SIAM  J.  Matrix  Anal.  AppL,  pp. 
844-881,  2008.  M.  W.  Mahoney  and  R  Drineas,  CUR  matrix  decompositions  for  improved  data  analysis.  Proc. 
Nat.  Acad.  ScL,  106:697-702,  2009.  S.  Wang  and  Z.  Zhang,  Improving  CUR  matrix  decomposition  and  the  Nystrom 
approximation  via  adaptive  sampling.  J.  Machine  Learning  Res.,  14:2729-2769,  2013. 
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kxn.  The  algorithm  is  novel  in  that  it  only  requires  one  pass  through  the  data  A  and  only  requires 
a  few  eolumns  of  A  to  be  in  memory  at  any  given  step  of  the  algorithm.  This  is  ideal  for  POD 
based  sehemes  sinee  the  high  dimensional  trajeetory  need  not  be  stored. 


To  illustrate  the  performanee  of  DEIM-CUR  as  eompared  to  a  leverage  seore  algorithm,  the  fol¬ 
lowing  example  builds  a  matrix  A  G  ]^300,000x300  j^^ving  the  form 


,  ^  1000 
7=1  ^ 


300  1 
T  ^  T 

*7^7  +  L  7*7^7’ 
7=11  J 


(3) 


where  G  and  y j  G  are  sparse  veetors  with  random  nonnegative  entries  (in  MAT- 

LAB,  Xj  =  sprand(30000, 1,0.025)  and  yj  =  sprand(300, 1,0.025)).  In  this  instantiation,  about 
18%  of  all  entries  of  A  are  nonzero.  The  form  (3)  is  not  a  singular  value  deeomposition,  sinee  the 
veetors  in  {x^}  and  {yj}  are  not  orthonormal;  however,  due  to  the  sparsity  of  these  veetors  this 
deeomposition  suggests  the  strueture  of  the  SVD:  the  singular  values  deeay  like  1  /  j,  and  with  the 
first  ten  singular  values  weighted  more  heavily  to  give  a  signifieant  drop  between  Oio  and  On. 

Figure  1 1  eompares  the  error  ||  A  —  CUR||  of  the  CUR  approximation  for  DEIM-CUR  to  two  meth¬ 
ods  based  on  leverage  seores.  For  the  first  method,  the  rank- A:  approximation  takes  the  k  rows  and 
eolumns  that  have  the  largest  leverage  seores  eomputed  using  only  the  leading  10  left  and  right 
singular  veetors;  the  seeond  method  is  the  same,  but  using  all  singular  veetors.  Both  methods  per¬ 
form  rather  more  poorly  than  DEIM-CUR,  whieh  elosely  traeks  the  optimal  value  O/t+i-  As  seen 
in  Figure  11,  the  DEIM-CUR  approaeh  delivers  an  exeellent  approximation,  while  seleeting  the 
rows  and  eolumns  with  the  leading  leverage  seores  does  not  perform  nearly  as  well. 


-^LS  (all)  — ^LS  (10)  — ^DEIM  —(Jk  +  i 


Figure  11:  Aeeuraey  of  CUR  approximations  using  k  eolumns  and  rows,  for  DEIM-CUR  and 
two  leverage  seore  strategies  for  the  sparse,  nonnegative  matrix  (3).  “LS(all)”  seleets  rows  and 
eolumns  having  the  highest  leverage  seores  eomputed  using  all  300  singular  veetors;  “LS(IO)” 
uses  the  leading  10  singular  veetors. 
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